perm filename LPC.F4[X,ALS] blob
sn#136799 filedate 1974-12-19 generic text, type C, neo UTF8
COMMENT ⊗ VALID 00006 PAGES
C REC PAGE DESCRIPTION
C00001 00001
C00002 00002 SUBROUTINE LPC(AIFFY,SPT,NPTS,M,NSP)
C00005 00003 SUBROUTINE INVFL(X,NP,M,RO,ERRN)
C00007 00004 SUBROUTINE FFTP(X,Y,M,L)
C00008 00005 SUBROUTINE RBITS(X,Y,M)
C00010 00006 SUBROUTINE USCRM(X,Y,M)
C00011 ENDMK
C⊗;
SUBROUTINE LPC(AIFFY,SPT,NPTS,M,NSP)
C THE PARAMETERS IN THE LIST ARE DEFINED AS FOLLOWS:
C AIFFY←THE INPUT DATA ARRAY OF REAL NUMBERS...A TIME SERIES
C NPTS←THE NUMBER OF POINTS IN THE TIME SERIES
C M← THE NUMBER OF FILTER COEFFICIENTS DESIRED.
C SPT← THE ARRAY IN WHICH THE REAL SPECTRUM IS RETURNED.
C NSP← THE NUMBER OF POINTS IN THE REAL SPECTRUM.
DIMENSION AIFFY(1),CF(50),SPT(1),X(512),Y(512)
TPI=2.*355./113.
10 NP=NPTS-1
C DON'T CALCULATE THE DIFFERENCE SIGNAL
DO 11 J=1,NP
11 X(J)=AIFFY(J)
C 11 X(J)=AIFFY(J+1)-AIFFY(J)
C APPLY A HAMMING WINDOW TO THE TIME SERIES.
DT=NP
X(NPTS)=X(NPTS-1)*.08
DO 12 J=1,NP
T=J-1
12 X(J)=X(J)*(.54-.46*COS(TPI*T/DT))
C CALCULATE AUTOCORRELATION COEFFICIENTS AND SOLVE
C FOR INVERSE FILTER COEFFICIENTS
100 DO 105 J=1,NPTS
105 Y(J)=X(J)
CALL INVFL(Y,NPTS,M,RO,ERRN)
C Y ARRAY RETURNS 1.,A1,A2,...,AM,0.,....0.
MM=M+1
DO 125 J=1,MM
125 CF(J)=Y(J)
300 DO 303 J=1,512
X(J)=0.
303 Y(J)=0.
MM=M+1
DO 310 J=1,MM
310 X(J)=CF(J)
C SHUFFLE DATA FOR CALC OF REAL TRANSFORM
MM2=MM/2 +1
DO 320 K=1,MM
K1=2*K-1
K2=2*K
X(K)=X(K1)
320 Y(K)=X(K2)
NMOD=1
323 IF (2**NMOD.GE.NSP) GO TO 326
NMOD=NMOD+1
GO TO 323
326 IF ((2*MM2).LT.MM) MM2=MM2+1
MOD=1
330 IF (2**MOD.GE.MM2) GO TO 340
MOD=MOD+1
GO TO 330
340 CALL FFTP(X,Y,NMOD,MOD)
CALL RBITS(X,Y,NMOD)
CALL USCRM(X,Y,NMOD)
C CALCULATE LOG OF RECIPROCAL INVERSE FILTER SPECTRUM
DO 350 J=1,NSP
350 X(J)=X(J)*X(J)+Y(J)*Y(J)
BIG=-1.0E+20
DO 360 J=1,NSP
IF (X(J).EQ.0) X(J)=1.*10E-36
SPT(J)=-10.*ALOG10(X(J))
360 IF(BIG.LT.SPT(J))BIG=SPT(J)
DO 377 J=1,NSP
377 SPT(J)=SPT(J)-BIG
RETURN
END
SUBROUTINE INVFL(X,NP,M,RO,ERRN)
DIMENSION X(512)
DIMENSION F(50),TF(50),A(50),R(50)
C CALCULATION OF M+1 LENGTH AUTOCORRELATION SEQUENCE
C FROM THE INVERSE FILTER INPUT SEQUENCE
MP1=M+1
DO 11 JJ=1,MP1
J=JJ-1
MMJ=NP-J
SS=0.
DO 10 I=1,MMJ
IPJ=I+J
10 SS=SS+X(I)*X(IPJ)
11 R(JJ)=SS
RO=R(1)
C FORTRAN IMPLEMENTATION OF LEVINSON'S METHOD FOR
C SOLVING THE AUTOCORRELATION MATRIX
F(1)=1.
ALPHA=R(1)
BETA=R(2)
A(1)=-R(2)/R(1)
GAMMA=A(1)*R(2)
DO 1 N=2,M
NM1=N-1
C=-BETA/ALPHA
IF (N-2) 2,2,3
3 DO 4 J=2,NM1
NN=N-J+1
4 TF(J)=F(J)+C*F(NN)
DO 5 J=2,NM1
5 F(J)=TF(J)
2 F(N)=C*F(1)
ALPHA=ALPHA+C*BETA
BETA=0.
DO 6 J=1,N
NN=N-J+2
6 BETA=BETA+F(J)*R(NN)
Q=-(R(N+1)+GAMMA)/ALPHA
DO 7 J=1,NM1
NN=N-J+1
7 A(J)=A(J)+Q*F(NN)
A(N)=Q*F(1)
GAMMA=0.
DO 8 J=1,N
NN=N-J+2
8 GAMMA=GAMMA+A(J)*R(NN)
1 CONTINUE
C CALCULATE NORMALIZED ERROR ERRN
C SM=0.
C DO 77 J=1,M
C 77 SM=SM+A(J)*R(J+1)
C ERRN=1.+SM/R(1)
C PLACE COEFFICIENTS IN PROPER LOCATIONS OF THE
C REAL VECTOR X FOR FFT-ING TO GET INVERSE SPECTRUM
DO 51 J=1,M
51 X(J+1)=A(J)
X(1)=1.
MP1=MP1+1
DO 123 J=MP1,512
123 X(J)=0.
RETURN
END
SUBROUTINE FFTP(X,Y,M,L)
DIMENSION X(512),Y(512)
N=2**M
L2=2**L
DO 1 LO=1,M
LMX=2**(M-LO)
LMM=LMX
LIX=2*LMX
SCL=6.283185/LIX
C TEST FOR PRUNING
IF (LO-M+L) 20,30,30
20 LMM=L2
30 DO 1 LM=1,LMM
ARG=(LM-1)*SCL
C=COS(ARG)
S=SIN(ARG)
DO 1 LI=LIX,N,LIX
J1=LI-LIX+LM
J2=J1+LMX
T1=X(J1)-X(J2)
T2=Y(J1)-Y(J2)
X(J1)=X(J1)+X(J2)
Y(J1)=Y(J1)+Y(J2)
X(J2)=C*T1+S*T2
1 Y(J2)=C*T2-S*T1
RETURN
END
SUBROUTINE RBITS(X,Y,M)
DIMENSION X(512),Y(512),L(9)
EQUIVALENCE(L9,L(1)),(L8,L(2)),(L7,L(3)),(L6,L(4)),
1 (L5,L(5)),(L4,L(6)),(L3,L(7)),(L2,L(8)),(L1,L(9))
C PERFORM BIT REVERSAL
DO 70 J=1,9
L(J)=1
IF (J-M) 71,71,70
71 L(J)=2**(M+1-J)
70 CONTINUE
JN=1
DO 60 J1=1,L1
DO 60 J2=J1,L2,L1
DO 60 J3=J2,L3,L2
DO 60 J4=J3,L4,L3
DO 60 J5=J4,L5,L4
DO 60 J6=J5,L6,L5
DO 60 J7=J6,L7,L6
DO 60 J8=J7,L8,L7
DO 60 JR=J8,L9,L8
IF (JN-JR) 61,61,62
61 R=X(JN)
X(JN)=X(JR)
X(JR)=R
FI=Y(JN)
Y(JN)=Y(JR)
Y(JR)=FI
62 JN=JN+1
60 CONTINUE
RETURN
END
SUBROUTINE USCRM(X,Y,M)
DIMENSION X(512),Y(512)
I=2**M
LO2=I/2
SA=3.141593/I
X(1)=2.*(X(1)+Y(1))
Y(1)=0.
LLL=LO2+1
X(LLL)=2.*X(LLL)
Y(LLL)=-2.*Y(LLL)
DO 33 K=2,LO2
J=I-K+2
ANG=SA*(K-1)
C=COS(ANG)
S=SIN(ANG)
AA=X(K)+X(J)
BB=Y(K)-Y(J)
CC=X(K)-X(J)
DD=Y(K)+Y(J)
XR=C*DD-S*CC
XI=S*DD+C*CC
X(K)=AA+XR
Y(K)=BB-XI
X(J)=AA-XR
Y(J)=-BB-XI
33 CONTINUE
RETURN
END